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ABSTRACT 

We report far-IR Herschel observations obtained between 70 pm and 500 /jm of two star-forming dusty condensations, Bl-bS and Bl-bN, in the 
Bl region of the Perseus star-forming cloud. In the Western part of the Perseus cloud, Bl-bS is the only source detected in all of the 6 PACS and 
SPIRE photometric bands without being visible in the Spitzer map at 24 /im. Bl-bN is clearly detected between 100 pm and 250 pm. We have 
fitted the spectral energy distributions of these sources to derive their physical properties, and find that a simple greybody model fails to reproduce 
the observed SEDs. At least a two-component model, consisting of a central source surrounded by a dusty envelope is required. The properties 
derived from the fit, however, suggest that the central source is not a Class object. We then conclude that while Bl-bS and Bl-bN appear to be 
more evolved than a pre-stellar core, the best-fit models suggest that their central objects are younger than a Class source. Hence, they may be 
good candidates to be examples of the first hydrostatic core phase. The projected distance between Bl-bS and Bl-bN is a few Jeans lengths. If 
their physical separation is close to this value, this pair would allow the mutual interactions between two forming stars at a very early stage of their 
evolution to be studied. 

Key words. Stars: protostars - Star: individual: [HKM99] Bl-bS - Star: individual: [HKM99] Bl-bN - Star: individual: [EYG2006] Bolo 81 - 
Star: individual: [EDJ2009] 295 



1. Introduction 

During the formation of a low-mass star, the first hydrostatic core 
(FHSC) phase is one of the shortest phases, lasti ng only 10 2 - 
10 3 yrs, depending on how the collapse proceeds ( Bate |20TT 



The FHSC phase is characterized by the presence of a central ob- 
ject with the mass of a giant planet and having a size of just a few 
astronomical units. This phase lasts until the central object tem- 
perature reaches ~2000 K. At that point, the molecular hydrogen 
dissociates, causing a second collapse that brings the central ob- 
ject to the typical size of a protostar, starting the Class phase 
( |Andre etaL 1993| l. The short duration of the FHSC phase makes 
it difficult to observe an FHSC in a star-forming region. 

Before the FHSC phase, the star-forming core (a "pre-stellar 
core") has a spectral energy distribution (SED) that can be mod- 
elled as a greybody spectrum. After this FHSC phase, the SED 
in the far-infrared (FIR) still resembles a greybody spectrum, but 
the emission of the central forming protostar starts to be visible 
at shorter wavelengths, i.e., A < 70 pm. Theoretical models (e.g., 
Omukai |2007||Saigo & Tomisaka|20TT) show that the SED de- 



viates from a pure greybody shape during the FHSC phase at 
wavelengths A ^ 200 /mi. At even shorter wavelengths, e.g., 
for A £ 30 pm, the emission is still too faint to be detected. 
Taking into account that the photometric bands ot the Spitzer 
MIPS instrument ( |Rieke et al.|2004[ > were centred in the FIR at 
24 pm, 70 pm, and 160 pm, the expected signature for a 1 M 



* Herschel is an ESA space observatory with science instruments 
provided by European-led Principal Investigator consortia and with im- 
portant participation from NASA. 



FHSC SED includes a detection at 70 pm and the non-detection 
at 24 pm (Commergo n et al. |2012| >. 

None of the FHSC candidates meets these photometric re- 
quirements. In the Perseus star forming region, the candidate 
Per 58 ( |Enoch et aLl[20T0l |Dunham etli^[20TT] l is visible 
in the 24 pm Spitzer map. Similarly, the candidate MMSI in 
Chamaeleon ( |Belloche et~al~||2Tj(36) |Cordiner et~aT7~|[20T2| > is 
observed both at 24 pm and at 70 pm. The proposed FHSC 
candidates L 1451 -mm ( |Pineda et aL"1|201 l|i, B olo 45 ( |Schnee| 



|et al.||20T2) , and CB 17 ( |Chen et al. ||2012| >, are indeed not 



visible at 24 pm, but they remain undetected also at 70 pm in 
the Spitzer maps, making the SED-based classification less firm. 
Finally, L1448-IRS2E ( |Chen et al. |2010| l was not even detected 
at 160 pm. It must be said, however, that the knowledge of the 
SED alone is not enough to firmly identify an FHSC. Spectral 
lines and interferometric observations are necessary and, indeed, 
these FHSC candidates were not proposed based on continuum 
data only. 

Further, the (non)detection in a certain band clearly depends 
on the distance to the source. To predict the expected flux from 
an FHSC, the distance of the star forming region in Taurus 



(150 pc) is often taken as reference distance (e.g., Tomida et 
al. |2010 Commer gon et al. |2012| l. Moving to larger distance, 



however, makes stronger the requirement on non detection at 
24 pm, while the detection at 70 pm may become less strin- 
gent. In fact, the non-detection at 70 pm has been explained by 
|Enoch eta l. (2010) on the basis that an FHSC could, in principle, 
be detected with Spitzer, but the "Cores to Disks" (c2d) survey 
( Evans et al.|2003[ ) was not sufficiently sensitive to detect these 
sources. The recent launch of the Herschel Space Observatory 



1 



Stefano Pezzuto et al.: Bl-bS and Bl-bN: two first hydrostatic cores candidates in the Perseus star-forming cloud 



satellite (Pilbr att et al.|2010j l has now opened up the possibility 
to observe in the FIR bands between 70 fim and 500 fim, with 
unprecedented sensitivity and spatial resolution. Among the ob- 
serving programs, the Herschel Gould Belt survey (GBS; Andre 



et al. 2010) aims to obtain a complete census of pre-stellar cores 



and Class sources in the closest star-forming regions. Since 
FHSCs can be considered extremely young Class sources, we 
expect to find new FHSC candidates among the numerous ob- 
jects discovered with Herschel. 

One of the targets clouds of the GBS is the star-forming 
region in Perseus molecular cloud, at an average distance 
of ~235 pc (Hirota et al. 2008)>. Perseus hosts low and 



intermediate-mass young stellar objects, and so it is roughly be- 
tween a low-mass star-forming cloud like Taurus, and a high- 
mass star-forming cloud like Orion. 

As a starting point for our analysis of these observations, we 
looked for sources detected in the Herschel map at 70 /mi with- 
out any counterpart in the corresponding Spitzer 24 /mi map. We 
found one such source and in this paper we discuss the physical 
properties that can be derived from the analysis of its SED. This 
object is spatially associated with the source Bl-bS discovered 
by Hirano et al. |(|1999]l. They found this source through interfer- 



ometric observations carried out with the Nobeyana Millimeter 
Array (NMA) at 3 mm. They also found a second source, Bl-bN, 
~20" north of Bl-bS, clearly visible in the Herschel PACS maps 
at A > 100 /mi. Combining the NMA data with other single dish 
observations bet ween 350 urn and 2 m m, and spectral observa- 
tions of H 13 CO + , [Hirano et al.~l ( |l999t concluded that these two 
sources were younger than already known Class objects. Their 
observations, however, were all taken in the Rayleigh-Jeans part 
of their SEDs, causing a large uncertainty in the determination 
of the temperature and, in turn, of the mass. Herschel data now 
extend the knowledge of the SEDs exactly where the peak of the 
emission falls. 

On the basis of Herschel fluxes, which are incompatible with 
a greybody, we reconsider the physical properties of these two 
sources, and we propose they are new promising FHSC candi- 
dates. The rest of the paper is so organized: in Section 2, the 
observations are presented and the algorithm of source detec- 
tion, CuTEx, is briefly described. In Section 3, we fit the SEDs 
with a two-component model that describes the FHSC emission 
in a simple way. We show that alternative approaches fail to 
reproduce the SEDs. In Section 4, we discuss the results, and, 
in Section 5, our work is summarized. Appendix A gives a de- 
tailed description of the SED fitting, including ancillary expla- 
nations on the colour corrections and on the Gaussian fitting 
of sources detected with Herschel. We also present a new for- 
mula to derive the mass of an object whose SED is fitted with a 
greybody. Appendix B describes an alternative approach to the 
sigma-clipping algorithm for outlier detection, where the thresh- 
old for the detection depends on the size of the sample. 

2. Observations and data analysis 

The observations of Perseus cover a total area of about 
10 deg 2 , mapped with the PACS (Photodetector Array Camera 
& Spectrometer; |Poglitsch et al.|2010|l and SPIRE ( Spectral and 
Photometric Imaging Receiver; (Griffin et al.|2010) |Swinyard et| 
|al.|2010| l instruments using the parallel mode with the telescope 
scanned at 60"/s. In this observing mode, the two instruments 
observe almost the same field with SPIRE delivering data in its 
three bands (250 /mi, 350 /mi, and 500 fim), and PACS in two out 
of three bands, 70 fim and 160 /mi in our case. A smaller area, 
6.7 deg 2 , was mapped again with PACS at 100 fim and 160 fim 



only, with the telescope scanned at 20"/s, to reach a higher spa- 
tial resolution and a better sensitivity. 

The observations we consider in this paper cover a region of 
about 5.6 deg 2 , centred on NGC 1333. Table [l] summarizes the 
observations. Preliminary results from the parallel mode obser- 
vations were recently presented by Sadavo y" et al.| ( 20T2) l, and 
by Bressert et al. (submitted). A complete description of this re- 
gion based on Herschel data will be given in an upcoming paper 
(Pezzuto et al., in preparation). 

Data were processed with version 8.0.1559 of HIPE ( |Ott| 
[| |2010| l, except: a) deglitching, which was performed with a 
sigma-clipping algorithm with a variable threshold detailed in 
Appendix B, and b) map reconstruction, which was made us- 



ing the ROMAGAL code (Trafica nteet al.|201 1| >. The zero point 
of the flux calibration for extended emission was established by 
comparing Herschel data with Planck and IRAS data over the 
same area ( |Bernard et al. |2010| . The result of the comparison is 
a set of offsets, one for each band, algebrically added to the ob- 
served images. For photometric measurements, however, where 
a background is estimated and subtracted from the source flux, 
as in this work, the zero points are actually not important. 

2. 1 . Source detection and photometry 

To identify sources in the field, we adopted the CuTEx algo- 
rithm ( |Molinari et al.|20TT i. This method computes the second 
derivative of the input image along four different directions. In 
the "curvature" image so obtained, we look for pixels where the 
second derivative is negative and its absolute value is above an 
user defined threshold with respect to the local fluctuactions. 
Isolated high-curvature pixels are rejected, and only groups of 
contiguous pixels, whose number depends on the sampling of 
the beam, are further considered. The curvature in each group is 
analyzed to verify the statistically significance of the local peak. 
The detection is done independently at each wavelength. Once 
a preliminary list of candidate sources is derived, photometry is 
recovered by fitting a 2D Gaussian profile. The position of the 
peaks of the candidate sources and their sizes, as derived from 
the curvature map, are used as initial guesses for the parameters 
of the 2D Gaussian. If two or more sources are closer than twice 
the width of the point spread function (PSF) at a given wave- 
length, their Gaussians are fitted simultaneously. 

For A > 160 fim, the diffuse emission contributes a strong, 
variable signal, across all spatial scales. An important issue is 
then how to estimate the contribution of the diffuse background 
at the source position. CuTEx models such emission by fitting a 
linear 2D-polynomial (i.e., a plane) inside a subregion centred on 
each source, of size 6 times the PSF at the specific wavelength. 
The plane and the Gaussian(s) are fitted simultaneously. There 
is, however, no general consensus on how the background can 
be estimated: different approaches can give different results and 
one of the main, and to large extent unknown, uncertainties in the 
derived fluxes, and in turn on the shape of the SED, comes from 
the background subtraction. For this reason is important to com- 
pare the results of flux extraction with more than one method. To 
this aim, we also used getsources (Men'shchiko v et al.|2012| l to 
derive the fluxes of our sources. 

Among the sources detected at 70 fim, we found that only 
one was not detected in the corresponding Spitzer 24 fim map. 
This source is one of a pair of very young objects first discovered 
at millimetre wavelengths by |Hirano et al. | ( | 1999[ > and dubbed 
Bl-bS and Bl-bN. Our 70 fim source corresponds positionally 
to Bl-bS, while Bl-bN is clearly detected between 100 fim and 
250 /mi. The separation between the two sources is about 20", 
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Table 1. The log of the obervations. Column OBSID reports the Observational Identifier which specifies an observation in the 
Herschel Science Archive. Each field was observed twice along two almost orthogonal directions to better remove the instrumental 
l/f noise. Date refers to the start of the observation; OD is the operational day (OD 1 is 14 May 2009). Centre gives the central 
coordinates of the map. Size is in arcminutes. The Observing mode column reports Parallel, for PACS and SPIRE data obtained 
in parallel mode, or PACS for data obtained with PACS only. PACS bands column reports the effective wavelength of the selected 
PACS filters. Speed is the scanning velocity of the telescope. 



OBSID 


Date 


OD 


Centre 


Size 


Obs. mode 


PACS bands 


Speed 


1342190326 
1342190327 


09/02/2010 
09/02/2010 


271 
271 


3 h 29 ra 39 s +30 d 54'32" 
3 h 29 m 41 s +30 d 53'58" 


135'xl50' 
150'xl35' 


Parallel 
Parallel 


70 pm, 160 fim 
70 fim, 160 fim 


60"/s 
60"/s 


1342227103 
1342227104 


22/08/2011 
23/08/2011 


831 
831 


3 h 30 ra 17 s +30 d 48'05" 
3 h 30 m 17 s +30 d 47'59" 


135'xl35' 
135'xl35' 


PACS 
PACS 


WOpm, 160yum 
100 pm, 160yum 


20"/s 
20"/s 



larger than the beam of PACS in all bands. It is also larger than 
the beam of SPIRE at 250 pm (~ 1 8"), but smaller than the beams 
at 350 pm and 500 pm (~25"and ~36", respectively). We at- 
tempted to fit two Gaussians to the SPIRE data using the posi- 
tions of Bl-bS and Bl-bN at 160 pm as initial guesses. During 
the fitting procedure the initial distance between the two objects 
was kept constant. 

Herschel data were complemented with archival Spitzer 
data from the c 2d survey (|Evans et al.||2003| l, SCUBA Legacy 
Catalogue data (Di Francesco et al. 2008 ) at 450 /im and 850 pm, 
and the Bolocam data ( |Enoch et al. 2006) at 1.1 mm. We ran 
CuTEx on all these maps except the Spitzer ones (see below). 
Again, the positions of the sources at 160 fim were used as initial 
guesses. The effective FWHM resolutions of the SCUBA Legacy 
Catalogue data at 450 fim is 17'.'3, so that the two sources can be 
considered to be fairly separated. At 850 fim, however, the effec- 
tive FWHM is 2279 and encompasses both Bl-bS and Bl-bN, 
causing the photometry to be more uncertain. At 1.1 mm, the 
Bolocam beam of 31 "makes it impossible to distinguish the two 
objects. 

Spitzer source S295 is very close to Bl-bS. It is visible in 
Herschel maps at A < 100 fim. S295 is commonly associated 
with the millimeter- wavelength source Bolo 81 ( |Enoch et al. 
2006). Since Bl-bS is close to the wings of the Spitzer PSF of 
S295, it could be that Bl-bS is not seen at 24 pm because of 
its proximity to the bright S295. To test this possibility, we used 
DAOPHOT on the Spitzer map to perform PSF photometry of 
S295, a method only possible due to the low background level at 
24 pm. Bl-bS remains undetected also after S295 is subtracted; 
from the analysis of the detections, we concluded that any source 
with a flux >0.2 mJy, corresponding to a >5 cr detection, should 
have been detected, and so we assumed 0.2 mJy flux as an upper 
limit at 24 pm for both Bl-bS/bN. At shorter wavelengths, the 
distance between the Bl-b sources and S295 is large enough and 
the background is so low and smooth, that to estimate an upper 
limit for the fluxes in the IRAC bands we just computed the rms 
of the background over a region of 11x11 pixels, centred at the 
160 pm positions. We used this value as the 1 <x upper limit for 
the peak flux which, along with the instrumental FWHM, can 
then be used to derive an upper limit at the IRAC wavelengths 
for a 2D Gaussian profile. The upper limits for Bl-bS are sys- 
tematically larger than those of Bl-bN. Since the former is close 
to the brighter S295, it is possible that the Bl-bS upper lim- 
its suffer some flux leakage from S295. For this work, however, 
shifting the upper limits by some percentage does not impact the 
fitting result (see below). The full set of measured fluxes is re- 
ported in Table [2] Image cut-outs at all wavelengths centred on 
Bl-bS, L5xl!5 in size, are shown in Fig.[T] 

The positions of the Bl-b sources from 24 pm to 1.1 mm are 
shown in Fig. [2] The triangle denotes the position of S295 taken 



Table 2. Photometric fluxes of B 1 -bS and B 1 -bN. The four rows 
denoted as IRAC/ report the 1 cr upper limits derived on the maps 
at the positions of the sources at 160 pm. Size is the instrumental 
FWHM. The row denoted as MIPS 1 reports the 5 cr upper limit 
at 24 pm, found by analyzing the Spitzer map with DAOPHOT 
inside an aperture of 12". For all the other wavelengths, the 
fluxes were measured by fitting a 2D Gaussian at the positions 
found by CuTEx. Size is the FWHM of the elliptical Gaussian 
fitted to each source, circularized and deconvolved with the in- 
strument beam size. The 70 pm flux of Bl-bN is equal to that 
of a point source having a 1 cr rms peak flux, so we consider it 
as an upper limit. The 160 pm data are from the map scanned 
at 20"/s. The uncertainties were estimated as the differences in 
the fluxes extracted with two independent methods: CuTEx and 
getsources (Men'shchikov et al.|2012 1. Herschel fluxes have had 
Gaussian fit corrections applied but no colour corrections, see 
Appendix A. 1 . 







Bl-bS 




Bl-bN 




A 


Instr. 


Flux 


Size 


Flux 


Size 


(pm) 




(Jy) 


(") 


(Jy) 


(") 


3.8 


IRAC1 


< 1.2 10-" 


1.66 


< 1.0 io- b 


1.66 


4.5 


IRAC2 


< 2.3 10~ 6 


1.72 


< 1.6 10" 6 


1.72 


5.8 


IRAC3 


< 9.4 10" 6 


1.88 


< 6.0 10~ 6 


1.88 


8.0 


IRAC4 


< 1.5 10 s 


1.98 


< 7.7 10~ 6 


1.98 


24 


MIPS1 


< 2.0 10~ 4 


12 


< 2.0 10~ 4 


12 


70" 


PACS 


0.22+0.14 


4.4 


<0.050 


6.9 


100" 


PACS 


2.29±0.26 


3.4 


0.61±0.31 


8.5 


160" 


PACS 


9.1 ±1.2 


5.8 


3.24±0.80 


8.9 


250 


SPIRE 


14.4 ±1.4 


10.4 


9.49±0.79 


14.3 


350 


SPIRE 


16.9 ±4.8 


17.8 


12.7 ±1.3 


18.4 


450* 


SCUBA 


19.1 ±4.6 


10.1 


13.5 ±3.7 


10.2 


500 


SPIRE 


15.8 ±7.8 


34.3 


14.6 ±4.9 


34.3 


850* 


SCUBA 


3.8 ±1.3 


19.8 


3.3 ±1.4 


22.3 


1100 r 


Bolocam 


3.67±0.26 Jy within 59'.'9x40'.'3 





Notes. Fluxes after the PSF correction discussed in Appendi x A; ^ Size 
refers to the inner part of the SCUBA beam, see Appendix A. 2; Enoch et al. 
(2006) : 2.61 Jy within 53"x59". 



from the Spitzer c2d catalog (Ev ans et al.| [2003), which agrees 
with our S295 positions at 70 pm and 100 pm. The positions 
of the Bl-b sources given in Hirano et al. |(|1999|l are pr ecessed 
to J2000. The position of Bolo 8 1 is from|Enoch et al.| ((2006] >. 
This source was later detected by Rosolowsky et al. '{ 2008 1 in a 



survey of ammonia cores, and dubbed NH3SRC 123. The cross 
labelled CuTEx Bolocam is the position we find using CuTEx on 
the 1 . 1 mm map. We also report in the figure the positions of the 
"core" and "outflow" recently found by Oberg et al. ( 2010 1 from 
CH3OH observations. The position they name "core" is halfway 
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between Bolo 81 and Bl-bN, while their "outflow" is slightly SE 
ofBl-bS. 

To derive the Herschel coordinates of the Bl-b sources we 
proceeded as follows. First, we averaged the PACS positions 
of S295 and compared them with the coordinates reported in 
the c2d catalog, finding an offset of Aa = 1'.'4 and A6 = 077. 
Assuming that the c2d coordinates are more accurate, these off- 
sets can be considered estimates of the absolute value of the 
PACS pointing errors in our maps. Next, we averaged the PACS 
coordinates of Bl-bS after adding the offset found for S295, 
finding a = 3 h 33 m 21'/3, 6 = +31 d 7'2774, with a total un- 
certainty (1 cr) of l'.'l. We did not average the SPIRE posi- 
tions because at 350 /zm and 500 /mi the two sources are more 
blended and the Gaussian fits are less reliable. We used the 
SPIRE 250 /zm position to estimate the offset between PACS 
and SPIRE: Aa = 472 and A6 = 370 Finally, for Bl-bN, we 
applied to the PACS 100 /mi and 160 /mi coordinates the same 
offsets found for S295; and to the SPIRE position at 250 /zm the 
offsets found between PACS and SPIRE for Bl-bS. The three 
pairs of coordinates were averaged finding a = 3 h 33 m 2172, 
6 = +31 d 7'4472, with a total uncertainty (1 cr) of 377. Note that 
all the offsets are much smaller than the separations between the 
objects, so that sources confusion due to a mis-pointing is ex- 
cluded. 

Even if it was already noted in the past that S295 and Bl-bS 
are spatially separated, in the literature there is some tendency to 
confuse the two sources. From Figs. [T] and [2] it is clear that these 
two objects are different. The association of Bolo 81 with S295 is 
very dubious since the latter is not detected longward of 100 /mi. 
The large offset, ~20", between Bolo 81 and the Bl-b sources 
also casts some doubts on the possible association between these 
objects. 



3. Results of SED fitting 

The simplest model that can be fitted to the data is an optically 
thin, isothermal greybody, which is appropriate for starless and 
pre-stellar cores 



31'07'5Cr 



k v B( T, v) 
D 2 



M 



where k v is the opacity (see the discussion in Appendix A. 3), 
B(T, v) is a blackbody at temperature T, D is the distance, M 
is the mass. Following the approach used for other GBS oberva- 



tions (e.g., Konyves et al. 2010[ ), the dust opacity index was fixed 
to 2. The best-fit models, considering only data at A > 160 /mi, 
have T = 9 K, M = 7.3 M Q for Bl-bS, and T = 8 K, M = 9.4 M Q 
for Bl-bN. 

We had to exclude from the greybody fitting the fluxes at 
short wavelengths because in a two-colour diagram of [100-160] 
vs. [250-350], the two Bl-b sources have colours that are bluer 
than those of a blackbody. Such a colour could only be repro- 
duced with a greybody given an unrealistic situation where the 
dust opacity increases with (far-IR/submm) wavelength. If the 



entire SEDs are fitted with a greybody, as in Equation (A.l i and 
leaving all the parameters free to vary, the best-fit models have 
B = 0, and T = 18 K for Bl-bS; and B = and T = 15 K, for 
Bl-bN. During the fitting procedure, we imposed the condition 
B > 0, and the value B — best matches the condition B < 0, 
corresponding to a dust opacity increasing with wavelength. 

Another way of fitting an SED with colours that are bluer 
than those of a blackbody is by assuming that the observed emis- 
sion is due to two components (see, e.g., Fig. 2 in |Pezzuto et] 



u 

Q 31"07'35 _ 
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Fig. 2. Positions (J2000) of the sources at all wavelengths; see 
text for references. Two black open circles are centred on the 
position of Bl-bS and Bl-bN after applying a positional off- 
set correction found by comparing the PACS and Spitzer posi- 
tions of S295, see text for details. The radius of the circles is the 
1 cr of the mean of the PACS coordinates, for Bl-bS, and of the 
mean of PACS 100 /zm, 160 /zm, and SPIRE 250 /zm coordinates 
for Bl-bN. The positions at 350 /zm, 500 /zm, and 850 /zm are 
very uncertain, because Bl-bS and Bl-bN are not resolved at 
those wavelengths. The blue open circle centred on Bolo 81 has 
a radius of 7", the pointing accuracy of Bolocam observations 
( |Enoch et al.|2006) . 

Table 3. The parameters of the best-fit models. T/, and Clb are the 
temperature and the solid angle of the blackbody; T g , Ao, B and 
Q. g are the parameters of the greybody (see Appendix A). Q/, and 
Cl g are given as radii ( y/£l/n). Uncertainties are 1 cr. The^ 2 are 
5.2 and 4.0 for Bl-bS and Bl-bN, respectively. For Bl-bN, the 
best-fit model has the smallest D./,, so we do not have an estimate 
for the lower uncertainty. 



Bl- 


T b 
(K) 


n b 

("AO) 


(K) 


*0 

Gum) 




Q g 
(") 


bS 
bN 


29.0 + " 
33.0+0 


c 9^+0.63 
J -0.50 
1 94+0.94 


8.0±0.5 
8.0±0.5 


100 -25 

120±45 


2.00±0.25 
2.00±0.25 


20.4+f 1 
17+"' 



|al. |20 02). For this reason, we fitted the Bl-b derived SEDs by 
adding the contributions of a blackbody embedded in a dusty en- 
velope whose emission is modelled with a greybody. We give a 
physical interpretation of this two-component model later on in 
this section. 

The details of the fitting procedure, as well as some addi- 
tional information on colour corrections applied to the PACS 
and SPIRE data, are given in Appendix A. The best result of 
the two-component fitting procedure is shown in Fig. |3]for both 
Bl-b sources, and the corresponding parameters of the best-fit 
models, Tb and Q^, for the blackbody component, and Tg, Ao, B, 
and for the greybody component, are listed in Table [3] 

The best-fit blackbody has a temperature of T ~ 30 K in 
both sources. Its size is about 130 AU at 235 pc for Bl-bS, and 
is much smaller in Bl-bN. For the latter, is not well con- 
strained since we do not have a 70 //m detection. The parameters 
of the greybody component, i.e., the dusty envelope, are similar 
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for both sources. The temperature is 8 K, not too much different 
from the 1 1 .7 K temperature found by Rosolowsky et al. ( 2008 ) 
for the ammonia core NH3SRC 123 detected at the position of 
B0I08I. 

From the relation m = acj/G, where c s is the sound speed, 
we can estimate the mass accretion rate. There are two limit- 
ing cases for the gravitational collapse of an isothermal sphere 
(see McKee & Ostriker 2007 , for a detailed discussion and refer- 
ences): in the fast collapse (the so-called Larson-Penston-Hunter 
solution) the infall is supersonic and a = 47, while in the slow 
collapse (or Shu's solution) the infall is sonic and a = 0.975. For 
T — 8 K we derive in the former case m = 5.3 x 10~ 5 M /yr, and 
m = 1.1 x 10~ 6 M /yr for the latter case. Taking into account that 
c s could include also contributions from a magnetic field or from 
turbulence, we can safely conclude that m < 1 x 10~ 4 M Q /yr. 

Turning these mass accretion rates into an accretion lumi- 
nosity is not immediate because we need to know the mass and 
the radius of the accreting source. The fact that we model the 
emission of the central object with a blackbody rules out the 
possibility to derive its mass. As far as the radius is concerned, 
the value we found from the fitting procedure appears larger than 
a few AU expected for an FHSC: 1 30 AU for B 1-bS and 48 AU, 
but with a 50% error, for Bl-bN. 

Such a large radius, however, can have a physical meaning. 



Bate (201 1 1 found that in a rotating envelope a disc with a size 
of ~ 100 AU can form before the stellar core; in the |Saigo &| 
Tomisaka ( 2011| Ts model, again for a rotating core, 100 AU 
correspond to the radius where the envelope's density profile 
changes from a r~ 2 power law to a flatten profile which ends on 
the core's surface, located at 20 AU when the FHSC first forms. 
The inner component of our model, then, can be physically inter- 
preted as describing the thicker and inner regions of the envelope 
that eventually merge with the compact central objects. In other 
words, the parameters Tb and Rh of the blackbody gives the av- 
erage properties of this hybrid component, without being able to 
distinguish where the envelope ends and where the core starts. 

The best-fit value for f3, the exponent of the dust opacity, is 
2, even if good fits for Bl-bN can be obtained also with /3 — 1.5; 
see Appendix A. The wavelength at which t = 1, Aq, is 100 /mi 
for Bl-bS and a value slightly larger for Bl-bN. 

The mass of the dusty envelope in Bl-bS is 9 ± 5 M , with 
the large uncertainty mainly due to the uncertainty in The 
mass of Bl-bN corresponding to the best-fit model is also 9 M Q , 
but the large uncertainty in Qj, and Q ? causes a large range of 
variation for the mass: 1 < M/M < 23. Previous estimates 



of the masses were made by Hirano et al. (T999J) who, with 
data at A > 350 /mi, derived for both sources smaller masses, 
M ~ 1.7 M , and higher temperatures, T ~ 18 K. IHatchell et al. 



(2007a) derived an upper limit for Bl-bS M 2 35 pc < 23.1 M Q , but 



using for the SED also the millimetre flux of Bolo 81 



The radius -y/Q g /^ of the envelope is ~20" for Bl-bS and 
~17" for Bl-bN. Such a radius, however, can not be immediately 
compared with the sizes given as FWHM in Table|2] which result 
from fitting the sources in the map with a 2D-Gaussian profile. In 
fact, when fitting the SED we assume that the source has a finite 
radius, R, while a Gaussian profile does not have any radius. To 
compare R with the observed FWHMs, we proceeded as follows: 
the integral 7 of a normalized 2D circular Gaussian over a circle 
of radius Rj is ( |Worz|2006| l 



Clearly, 1=1 only for R; — > 00, but already when R; = 3cr 
I ~ 0.993 ^ e men define the radius Rg of a 2D Gaussian 
as Re = 3<x, or, using the relation between FWHM and cr, 
Rg = 1.29xFWHM. For Bl-bS, the best-fit radius R is 2a.'4, 
a circular 2D-Gaussian with Rg — R has FWHM = 15'.'8, a value 
in between the measured size at 250 /mi and at 350 /mi. For Bl- 
bN, the best-fit radius is 17", which corresponds to FWHM = 
13", very close to the measured size at 250 /mi. 

It is evident from Table [2] however, that the size derived 
from the observations is not well defined, since the FWHM in- 
creases with wavelength. This trend is to some extent related to 
the instrumental beam size. For instance, the FWHM of Bl-bS 
is smaller at 100 //m than at 70 /mi, which is likely related to 
the smaller PSF of the intrument at 100 /mi (scanning speed of 
20"/s, nominal PACS compression mode) than at 70 /mi (scan- 
ning speed of 60"/s, double PACS compression mode). This 
trend is also likely to be related to the increase with wavelength 
of the background level, which makes it more difficult to disen- 
tangle the genuine source emission from the extended emission. 
It is physically plausible, however, to assume that this trend is, 
at least in part, intrinsic to the sources, given the fact that this 
phenomenon has been already observed elsewhere. 

In high-mass star-forming regions, for instance, an increase 
of the size of the dense cores with wavelength is known and is 
taken into account by scaling the flux proportionally to the ra- 
tio of measured radii, assuming the size at 160 /mi as the fidu- 
cial radius (see, e.g., Motte et al. 2010| Giannini et al.|20T2] l. 
This approach is justified by theoretical models showing that the 
emission can be described as coming from an almost isothermal 
envelope whose mass increases linearly with the radius. 

For cores in nearby star-forming regions, where the spatial 
resolution is high enough to resolve the structure of the enve- 
lope, the wavelength-radius relationship could be a consequence 
of temperature stratification in the envelope. Namely, at shorter 
wavelengths, we see the inner and warmer part of the envelope. 
At larger wavelengths, however, the emission comes from the 
outer and generally colder part of the core. In our model, part 
of the blackbody radiation is absorbed by the envelope which 
should cause a temperature stratification. We did not however 
treat the radiative transfer of the system, assuming instead, for 
simplicity, an isothermal envelope. The temperature T g so de- 
rived, and in general the whole set of parameters T g ,f3, Aq, Q g , 



along with the mass M from Equation ( |A.9[ ), can then be con- 
sidered to be the averaged properties of the outer regions of the 
envelope, in the same way as T/, and R/, describe average prop- 
erties of the inner regions as discussed above. In other words, 
while our objects likely have a temperature stratification which 
causes an inrease of the observed radius with wavelength, our 
model shows that the SEDs can be described in a simple way as 
having one radius and one temperature at all wavelengths. 

The consistency of our approach can be seen from the fact 
that: 1) the derived size is comparable to the observed sizes at 
A ~ 300 //m, and 2) Aq ^ 100 /mi. In fact, in a greybody model, 
the dependence on the solid angle disappears in the optically thin 
regime, so that the derived Q g values are most likely to represent 
the sizes corresponding to the largest wavelength where the en- 
velope is not yet completely thin. The robustness of our result 
can be verified by noting that T g and M are very similar to the 
values found with the isothermal greybody model discussed at 
the beginning of this section. 



7=1- exp ?<r 



1 More precisely, / = 0.9889..., slightly smaller than 0.9973, the in- 
tegral over the interval ±3<x for a ID-Gaussian. 
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Table 4. The luminosities of our sources, all in L . Lb is the lu- 
minosity of the blackbody, which in our model corresponds to 
the so-called internal luminosity; L em is the amount of black- 
body radiation that escapes the envelope, i.e., f B v (Tb)e~ Tv dv; 
L g is the luminosity of the greybody component; Lbd is the sum 
of L em and L g , the total observable luminosity predicted by the 
model; Lsed is the observed luminosity found by integrating the 
measured fluxes. 



Bl- 


u, 






Lboi 


i-SED 


bS 
bN 


0.45 
0.11 


0.26 
0.04 


0.22 
0.22 


0.48 
0.26 


0.49 
0.28 



Having modelled the observed SED with two components 
allows us to derive separately the internal luminosity L mt of the 
sources and the bolometric luminosity Lboi- 

In Table |4] we report: Lb, the luminosity of the black- 
body component. This is found with the standard equation 
Lb = (the parameters are those of Table pi with 

Rb = ( y/Clb/n)). At any frequency v, the envelope absorbs a 
part of the blackbody emission leaving the amount B v (Tb)e~ Tv 
free to escape from the envelope. The integral over frequency 
of this radiation is reported in the table as L em . Clearly, the 
fraction (1 - e~ Tv ) is absorbed inside the envelope. The column 
labelled L g reports the luminosity of the greybody component 
found by numerical integration of Equation (A.l i in the range 
1 /mi < A < 10 mm. The observable luminosity predicted by the 
model is then L em + L g , which we report in the column Lboi- It 
can be compared with the observed luminosity that we report in 
the column Lsed, found by integrating the measured fluxes. 

If there is no external source of energy, then by the conser- 
vation of energy L g — Lb - L em . The contribution of the inter- 
stellar radiation field to the sources luminosity could be, then, 
estimated as Lisrf = L g - (Lb - L em ). For the two sources, how- 
ever, we find quite different results: Lisrf = 0.03 L Q for Bl-bS 
and 0.15 L Q for Bl-bN. While a small difference between Lisrf 
in the two sources could be reasonable, such a large discrepancy 
is more likely due to the uncertanties in the best-fit parameters, 
in turn due to the uncertanties in the measured fluxes, or to the 
consequence of the non-proper treatment of the radiation trans- 
fer in our model. 

Finally, as an example of a model more rigorously treating 
the radiative transfer, we fitted the SEDs using the on-line fit- 
ting tool by Robitaille et al. (2007), that compares the observed 
fluxes with a grid of theoretical SEDs. This comparison is im- 
portant also because in this way we can test the hypothesis that 
our sources are actually more evolved than an FHSC. The mod- 
els by Robitaille et al. ( 2007 ) are indeed appropriate to describe 
the emission from an evolved young stellar object (YSO), where 
a central star already formed. On the contrary, their grid does not 
contain models that are appropriate for very young sources like 
an FHSC. We found that no model can account for both the up- 
per limits in the Spitzer bands and the far-IR/submm fluxes. The 
SEDs of our sources are not compatible with any of the 200 000 

( |2QQ7> . 



models of YSOs in the grid of Robitaille et al 



4. Discussion 

The following features can be extracted from the observed SEDs 
of our sources, prior to any modelling: the colours in the PACS 
bands are not compatible with those of a greybody, while the 
SEDs of the starless or pre-stellar cores can be modelled with 



Table 5. A comparison of the observed SEDs between 24 /mi 
and 160 /mi for the proposed FHSC candidates. MMSI is in 
Chamaeleon, CB 17 is a dark globule in the constellation of 
Camelopardalis, and all the other sources are in the Perseus star- 
forming region. PACS fluxes and upper limits are from this work 
w ith the exception of the PACS 100 /mi flux of CB 17, quoted 
in |Chen et al. | ( [20T2] l. Spitzer data are from the given reference 
with the exception of the 24 /mi upper limit for Bl-bS and Bl- 
bN, derived by us. For Bolo 45, Schne e et al.| ( |2012"} report that 
this source was not detected in the Spitzer map, but no upper 
limit is given. 
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a greybody. An additional blackbody is required, which hints at 
the presence of an inner and warmer compact component. This 
second component is not visible at A & 70 /mi, so it is not as 
evolved as other nearby protostellar objects. A further observa- 
tional signature, although indirect, of the very young age of Bl- 
bS comes from the statistics of the YSOs detected in Perseus. 
The preliminary analysis with CuTEx of the Western region of 
Perseus gives a list of ~ 40 tentative sources detected in all of the 
six Herschel bands. All of them are visible in the 24 /mi Spitzer 
map, the only exception being Bl-bS. This fact alone points to 
some peculiarity with this source. 

In Table [5] we report for all the FHSC candidates the photo- 
metric data available in literature, or derived by us, in the funda- 
mental range between 24 /mi and 160 /mi. As already stated in 
Section 1, the SED only is not enough to identify an FHSC. The 
other sources reported in Table [5] have been proposed as FHSC 
candidates based on other observational properties, e.g., the u-v 
visibilities from interferometric observations. On the other hand 
it is evident that as far as the SED is concerned Bl-bS remains 
an exception also when compared with the other candidates. It 
is the only source that shows a clear evidence of an SED more 
evolved than a pre-stellar core, but less evolved than a Class 
source. 

Other age indicators directly derived from the SEDs are the 
ratio L35o/Lb i, i.e., the luminosity at A > 350 /mi over the 
bolometric luminosity, and Lbd, i.e., the temperature of a black- 
body having the same mean frequency as the observed SED. 
For a Class source, L35 /L b oi > 5% ( [Andre et aT1|1993| l and 
Lboi < 70 K ( |Chen et al. |1995| l. Our sources have values quite 
extreme with respect to these limits, thus confirming that they are 
indeed very young. Bl-bS has L^q/L^x = 25% and Lb i = 18 K; 
Bl-bN has L35o/Lb i = 35% and L D oi = 14 K. From these sig- 
natures alone, however, it is not possible to distinguish an FHSC 
from a pre-stellar core, or a starless core. 

In summary, Bl-bS shows all the features expected in the 
SED of an FHSC. Bl-bN, the companion 20" north of Bl-bS, is 
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not detected at 70 fim, but otherwise it shares all of the charac- 
teristics of B 1 -bS . 

For both sources, only the bolometric luminosities, found by 
integrating the observed fluxes, do not fit with FHSC predic- 
tions: L bo i is 0.49 L and 0.28 L for Bl-bS and Bl-bN, re- 
spectively, in excess of the maximum foreseen for an FHSC, 
L $ 0. 1 L (|Omukai 2007 1. Such a high luminosity could be as- 



cribed to imperfect background removal, but even if we limit the 
integration to the PACS bands, where the diffuse emission is less 
prominent, we find for Bl-bS Lpacs = 0.14 L . It then seems 
reasonable to expect L ^ 0.2 L , a value that also excludes the 
possibility that Bl-bS is a very low luminosity objects (VeLLOs; 
e.g., |Di Francesco et al.|2007" i. 

The observed bolometric luminosity, however, can be larger 
than the internal luminosity due to external heating of the en- 
velope by the interstellar radiation field. We have seen in the 
previous section that is not easy to estimate Lisrf from our data. 
Instead, we can make use of the relationship between L mt and 
F-iq found by punhamet al. ( 2008| >. For B 1 -bS, their relationship 
gives Li nt ~ 0.04 L , a value which seems difficult to reconcile 
with the observed luminosities. On the other hand, Dun ham et| 
al. (2008 ) derived the relation by modelling the protostar with 
a"nxed temperature of 3 000 K and a luminosity in the range of 
0.03 — 10 L , which implies a radius in the range 0.5 — 11 R Q . 
Such values are reasonable for a Class source, but could be in- 



appropriate for a still-younger object, like a FHSC. [Dunham et 
al. ( 2008| Ts result, however, clearly implies that for our sources 
Lint is surely smaller than the measured L DO i. The precise factor, 
however, remains unknown. 

Recently, [Commergo n et aL~] (|2012) found that the luminos- 
ity of an FHSC becomes larger than 0.1 L only at late times, so 
that the vast majority of FHSCs appear as VeLLOs for most of 
all their lifetime. Nonetheless, they also found cases where only 
a few hundred years after the formation of the FHSC, the pre- 
dicted luminosity is as high as 0.3 L Q , similar to what we derive 
for Bl-bN. Such an increase in luminosity is due to the increase 
of mass of the inner core due to the accretion from the envelope. 

Another possible contribution to the luminosities of our 
sources is through contamination of the SEDs by an outflow, as 
recently found by |Maury~ et al. ( 20 1 2) for the prototypical Class 
protostar VLA1623. In our cases, however, the existence of an 
outflow in Bl-bS/bN is not completely clear, so it is not possible 
to explore the reliability of this hypothesis. Drabek et al. (2012) 
found that 12 CO line emission can contribute to the dust contin- 
uum at 850 fim. They estimated, however, that in NGC 1333, at 
positions far from outflows, the contamination is less than 20%. 
For the Bl region, Sadavoy et al. (in preparation) found that CO 
(3-2) emission appears to be relatively minor towards Bl-bS and 
Bl-bN, therefore any contamination of the dust continuum from 
the gas is likely insignificant at 850 //m. 

The results of the fits we performed must be interpreted cau- 
tiously given the uncertanties in the fluxes. We tentatively sug- 
gest, however, that Bl-bN is still slightly younger than Bl-bS. 
It has an envelope slightly more optically thick than that of Bl- 
bS, explaining the lack of detection in the shortest PACS band. 
For both sources, the envelope is cold and massive, and its in- 
ternal thermal pressure alone insufficiently high to be able to 
prevent its gravitational collapse. This conclusion comes from 
comparing the mass of the envelope, and the corresponding crit- 
ical Bonnor-Ebert mass m^E, the largest mass that an isothermal 
sphere of gas bounded by pressure can have without collapsing. 
When M env > m BE , the internal thermal pressure may not be 
high enough to support the core against internal gravity and ex- 
ternal pressure, so that the envelope is undergoing, or is about to 



collapse. We found for m BE , see Equation ( |A.3[ ), 0.4 M (Bl-bS) 
and 0.3 M ( Bl-bN). In both cases, M env » »ibe so that the en- 
velopes cannot be supported by internal thermal pressure alone, 
even if we cannot exclude that other physical mechanisms like 
the amount of turbulence or the strength of the magnetic field 
(see, e.g., Bas u et al.|2009[ ), may still contribute to the support 
of the cloud against the gravitational collapse 

As far as the magnetic field is concerned, however, the 
largest mass that a magnetic field of 31 //G (as estimated for 
the B 1 region by Matth ews & Wilson|20 02 ) can support is only 
0.11 M (from Stut z et al.|2007|) , assuming a radius of 0.023 pc 



(20"at 235 pc, see Sl g in Table |3|. 

To estimate the stability against turbulence we computed the 
virial parameter a (Bert oldi & McKee |1992| ) and we consider a 
core gravitationally bound if 

5cr 2 R 

a — s> 2 

GM 

cr is the velocity dispersion found by adding in quadrature 
the non-thermal component of the NH3 lines measured by 
Rosolowsky et al.| p508] ) for NH3SRC 123, and the thermal 
component of a mean particle of molecular weight /i = 2.33: 
o- NT = 0.325 km/s and cr T = 0.204 km/s. From Table [3]/? = 
7.19 x 10 16 cm and R = 5.99 x 10 16 cm for Bl-bS and Bl-bN, 
respectively, for a distance of 235 pc. Then, a = 0.44 and 0.37 
for the two sources. Clearly, the exact values of a are not well 
constrained given the large uncertainties, especially in the mass 
of the sources. As long as, however, M ^ 2 M G , a £ 2 

Important information on the evolutionary stage of a source 
can be obtained from the observations of its outflow. The con- 
clusions drawn from spectra l observations (e.g., IHatchell et al. 



2007b Oberg et al. 2010), however, have been affected by 



the lack of an adequate knowledge of the spatial distribution 
of the sources. Recall that Hirano et al. (1999) detected only 



Bl-bS/bN, while in the Bolocam survey at 1.1 mm ( |Enoch| 



|et al.||2006| only Bolo 81 was detected. Also, recall that the 
Spitzer maps ( |Rebull et al.|2007"] ) show only S295, while in the 
SCUBA map at 850 /jm Bl-bS and Bl-bN are blended. Thanks 
to Herschel data, for the first time we can see clearly Bl-bS/bN 
and S295 in the same map, in particular in the PACS bands. This 
new definition will allow a better understanding of the spectro- 
scopic observations. For instance, by comp aring our Fig.|2]with 
the CH3OH maps by |0berg et al. ( 2010| ), it appears possibile 



that Bl-bS coincides with the peak in the CH3OH map, in be- 
tween the two identified outflows (the authors give the coordi- 
nates of only one outflow, which we reported in Fig. [2}. 

The presence of outflows makes it possible that one of the 
two sources is actually a local density enhancement (knot). 
Indeed, [Gueth et aL~| ( |2003| ) concluded that emission knots gen- 
erated by the shocked outflow in the vicinity of a protostellar 
object can be misinterpreted as, e.g., starless clumps. This con- 
clusion, however, has been derived from mm and submm maps; 
to what extent it holds also in the FIR is unknown. 

We finally note that the number of FHSC candidates found 
in Perseus is quite large. |Pineda et al.| ( |201 \\ derived that the 
number of expected FHSCs in this region should be <0.2, 30 
times smaller than the 6 proposed candidates reported in Table[5] 
Beside the obvious argument that eventually none of the 6 can- 
didates may be confirmed to be an FHSC, there are other ways 
to explain this discrepancy. The expected number of FHSCs is 
found by assuming that the ratio of the number of FHSCs over 
the number of Class sources is equal to the ratio between their 
respective lifetimes. Since the two ratios are equal only in the 
case of continuous star formation, a first hypothesis is that we 
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are observing a burst of star formation in Perseus. The other pos- 
sibility, which like the previous one has already been suggested 
by Pineda et al.|(|201 1), is that we assume incorrect lifetimes for 



the FSHC phase, the Class phase, or both. Also the number 
of Class sources may be underestimated (Schnee e t al.|2012) , 
which would also impact the estimated lifetime. Though the cat- 
alog of detected sources in Perseus is not yet released (Pezzuto 
et al., in preparation), it is reasonable to suppose that Herschel 
observations will update the number of starless cores, pre-stellar 
cores and Class sources for all nearby star-forming regions. 
With this information, it is possible that the expected number of 
FHSCs will change. 

Finally, we note that our sources are quite close, with a pro- 
jected separation of 20". If the physical distance is not very 
different from the projected one, then Bl-bS and Bl-bN are 
^4700 AU apart, corresponding to few Jeans lengths at 10 K. 
It is then possible that these sources formed more or less at the 
same time from the fragmentation of a larger structure. This pos- 
sibility, if confirmed, would explain why we found two FHSC 
candidates at close proximity. 
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In this paper, we have presented the results of fitting the SEDs 
of two sources in the Perseus star-forming region. Data were de- 
rived from Herschel photometry observations collected within 
the Gould Belt Survey program (Andre et al. 2010 1. The two 
sources are Bl-bS and Bl-bN, discovered by |Hirano et al. 



( |1999| l from interferometric observations at millimeter wave- 
lengths. Bl-bS is the only source in the Western part of the star 
forming region in Perseus, detected in all six Herschel bands, 
and not visible in the Spitzer 24 /mi map. These two criteria are 
important for the SED-based detection of FHSC candidates. The 
SED alone is not enough to determine the evolutionary status of 
a source, and other proposed FHSC candidates have been ob- 
served with interferometry to better assess their status. But lim- 
iting to the photometric criteria Bl-bS is the only candidate that 
satisfies the detection at 70 pm and the non detection at 24 pm. 

Herschel data were complemented with Spitzer, JCMT- 
SCUBA, and CSO-Bolocam data. The resulting SED was fit- 
ted with a two-component model that describes the emission in 
terms of a blackbody, which roughly corresponds to the compact 
central object, and a greybody. i.e., the surrounding dusty enve- 
lope. We also tried to model the SED with a simple greybody 
alone, which is adequate for a pre-stellar core, and with a proper 
radiative transfert model (Robitaille et al. 2007 ), suited for more 
evolved sources. Both latter models failed to reproduce the SED 
over the whole observed spectral range. We conclude that Bl- 
bS shows almost all the expected characteristics expected in an 
FHSC. Only its luminosity is too high with respect to the theo- 
retical predictions. Indeed, such a high luminosity is at present 
the strongest argument against our candidate being an FHSC. 

The SED of Bl-bN is similar to that of Bl-bS, with the im- 
portant difference that the former is not detected at 70 pm. For 
the rest, this source seems similar to Bl-bS. The two sources are 
few 10 3 AU apart, corresponding to few Jeans lengths. It is then 
possible that these two sources formed at almost the same time 
from the fragmentation of a larger structure. 

Further observations at long wavelengths with high spatial 
resolution, such those that ALMA will undertake, are clearly 
needed to understand better the nature of B 1 -bS and B 1 -bN. It is 
reasonable, however, to conclude that Herschel has observed two 
young objects which could be in a very early stage of protostellar 
formation, and we propose them as two new FHSC candidates. 



Andrae, R. 2010, astro-ph:1009.2755v3 

Andre, P., Ward-Thompson, D., & Barsony, M. 1993, ApJ, 406, 122 
Andre, P., Men'shchikov, A., Bontemps, S. et al. 2010, A&A, 518, 102A 
Arzoumanian, D., Andre, P., Didelon, P. et al. 2011, A&A, 529, 6A 
Basu, S., Ciolek, G.E., Dapp, W.B. & Wurster, J. 2009, New A, 14, 483 
Bate, M.R., 2011, MNRAS, 417, 2036 

Beckwith, S.V.W., Sargent, A.I., Chini, R.S., & Gusten, R. 1990, AJ, 99, 924 

Belloche, A., Parise, B., van der Tak, EES. et al. 2006, A&A, 454, L51 

Bernard, J.-R, Paradis, D., Marshall, D.J. et al. 2010, A&A, 518, L88+ 

Bertoldi, F, & McKee, C.F. 1992, ApJ, 395, 140 

Chen, H., Myers, P.C., Ladd, E.F., Wood, D.O.S. 1995, ApJ, 445, 377 

Chen, X., Arce, H.G.., Zhang, Q. et al. 2010, ApJ, 715, 1344 

Chen, X., Arce, H.G.., Dunham, M.M. et al. 2012, ApJ, 751, 89 

Commercon, B., Launhardt, R., Dullemond, C, & Henning, Th. 2012, A&A, in 

press (astro-ph: 1207.0656) 
Cordiner, M.A., Charnley, S.B., Wirstrom, E.S., & Smith, R.G., 2012, ApJ, 744, 

131 

Di Francesco, J., Evans, N.J. II, Caselli, P. et al., 2007, in Protostars and Planets 
V, B. Reipurth, D. Jewitt, & K. Keil (eds.), University of Arizona Press, 
Tucson, 17 

Di Francesco, J., Johnstone, D., Kirk, H., MacKenzie, T. & Ledwosinska E. 

2008, ApJS, 175, 277 
Drabek, E., Hatchell, J., Friberg, P. et al. 2012, MNRAS, in press 
Dunham, M.M., Crapsi, A., Evans, N.J. II, et al., 2008, ApJS, 179, 249 
Dunham, M.M., Chen, X., Arce, H.G., et al., 201 1, ApJ, 742, 1 
Enoch, M.L.. Young, K.E., Glenn, J. et al. 2006, ApJ, 638, 293 
Enoch, M.L., Lee, J.-E., Harvey, P., Dunham, M.M. & Schnee, S. 2010, ApJ, 

722, L33 

Evans, N.J. II, Allen, L.E., Blake, G.A. et al. 2003, PASP, 1 15, 965 
Giannini, X, Elia, D., Lorenzetti, D. et al. , 2012, A&A, 539, 156A 
Griffin, M.J., Abergel, A., Abreu, A. et al. 2010, A&A, 518, L3 
Gueth, E, Bachiller, R., & Tafalla, M. 2003, A&A, 401, L5 
Hatchell, J., Fuller, G.A., Richer, J.S., Harries, T.J. & Ladd, E.F. 2007a, A&A, 
468, 1009 

Hatchell, J., Fuller, G.A. & Richer, J.S. 2007b, A&A, 472, 187 
Hildebrand, R.H. 1983, QJRAS, 24, 267 

Hirano, N., Kamazaki, T, Mikami, H., Ohashi, N, Umemoto, T. 1999, in Proc. 

of Star Formation 1999, ed. T. Nakamoto, Nobeyama Radio Obs., 181 
Hirota, T, Bushimata, T., Choi, Y.K. et al. 2008, PASJ, 60. 37 
Konyves, V., Andre, Ph., Men'shchikov, A. et al. 2010, A&A, 518, L106 
Matthews, B.C. & Wilson, CD. 2002, ApJ, 574, 822 
Maury, A., Ohashi, N., & Andre, P. 2012, A&A, 539, 130M 
McKee, C.F., & Ostriker, E.C. 2007, ARA&A, 45, 565 
Men'shchikov, A., Andre, P., Didelon, P. et al. 2012, A&A, 542, A81 
Molinari, S., Faustini, F, Schisano, E., Pestalozzi, M., & Di Giorgio, A. 2011, 

A&A, 530, 133A 
Motte, F., Zavagno, A., Bontemps, S., et al. 2010, A&A, 518, 77M 
Oberg, K.I., Bottinelli, S., J0rgensen, J.K. & van Dishoeck, E.F. 2010, ApJ, 716, 

825 

Omukai, K. 2007, PASJ, 59, 589 



8 



Stefano Pezzuto et al.: Bl-bS and Bl-bN: two first hydrostatic cores candidates in the Perseus star-forming cloud 



Ott, S. 2010, ASP Conference Series, 434, 139 

Pezzuto, S., Grillo, R, Benedettini, M. et al. 2002, MNRAS, 330, 1034 
Pilbratt, G.L., Riedinger, J.R., Passvogel, T. et al. 2010, A&A, 518, 1A 
Pineda, J.E., Arce, H.G., Schnee, S. et al. 201 1, ApJ, 743, 201 
Poglitsch, A., Waelkens, C, Geis, N. et al. 2010, A&A, 518, L2 
Rebull. L.M., Stapelfeldt, K.R., Evans, N.J. II et al. 2007, ApJS, 171, 447 
Rieke, G.H., Young, E.T., Engelbracht, C.W. et al. 2004, ApJS, 154, 25 
Robitaille, T.P., Whitney, B.A., Indebetouw, R. & Wood, K. 2007, ApJS, 169, 
328 

Rosolowsky, E.W., Pineda, J.E., Foster, J.B et al. 2008, ApJS, 175, 509 
Sadavoy, S.I., Di Francesco, J., Andre, P. et al. 2012, A&A, 540, A10 
Saigo, K. & Tomisaka, K. 201 1, ApJ, 728, 78 
Schnee, S., Di Francesco, J., Enoch, M. et al. 2012, ApJ, 745, 18 
Stutz, A.M., Bieging, J.H., Rieke, G.H. et al. 2007, ApJ, 665, 466 
Swinyard, B.M., Ade, P., Baluteau, J.-P. et al. 2010, A&A, 518, L4 
Tomida, K., Machida, M.N., Saigo, K., Tomisaka, K., & Matsumoto, T. 2010, 
ApJ, 725, L239 

Traficante, A., Calzoletti, L., Veneziani, M. et al. 201 1, MNRAS, 416, 2932 
Worz, S. 2006, PhD thesis, University of Hamburg. Aka, Berlin, Germany 



Once the minimum^ 2 was found, the 1 cr uncertainties in the 
parameters were found by consi dering all the m odels with a% 2 in 
the range x 2 mia < X 2 ^ X 2 min + 1 ( jAndrae |2010| . For Bl-bS, only 
8 models were found in this x 2 range. For Bl-bN, not having a 
70 fim flux had the consequence that 242 models, out of a total of 
27646, had y 2 < y 2 +1. Among these 242 models, the best fit 
was the one reported in Table[3] but with very large uncertainties. 
To decrease the uncertainties, we added another constraint and 
selected only the 107 models with Aq < 200 fim, in analogy with 
the results of the selection for Bl-bS. Finally, assuming that the 
dust properties are the same in the two sources, we considered 
only the models with B — 2, discarding 70 models with B — 1.5. 

For the Bonnor-Ebert mass, we used the formula reported in 
Konyves et al^pJTO) 

Ra 2 

m BE = 2.4— (A. 3) 



Appendix A: SED fitting and photometry 

For our two-component models, we modelled each set of data 
as the sum of a blackbody and a greybody envelope: F v = 
B{T,,)e- T Cl b + G(T g ,A ,B)Q. g , where 



G(T g ,A ,B) = (1 - e- T )B v (T g )Cl g 



(A.l) 



with t, the optical depth, parametrized as a power-law, i.e., 



(A.2) 



where Aq is the wavelength at which t = 1, and £l g and Q/, are 
the solid angles of each respective component. The independent, 
unknown parameters are then Tb,T g ,Ao, and B while the solid 
angles are computed by scaling the model fluxes to the data. 
This model describes our sources in terms of two components: 
a central source, the blackbody, embedded in a dusty envelope 
whose emission is modelled with a greybody. Since we do not 
impose that the envelope is optically thin at all wavelengths, the 
blackbody radiation is attenuated by a factor that describes the 
absorption due to the dust opacity. The absorption term, strictly 
speaking, implies that the envelope cannot be isothermal, so, as 
explained in Section 3, the set of parameters T g , Aq,B, Q. g should 
be considered only as describing the average physical conditions 
of the envelope. The density profile in the cores usually follows 
a power law, p cc r~ q , with typically q S> 2, so that the mass, 
and then the flux, increases with the radius. T g should be then 
indicative of the temperature in the outer region of the envelope. 

The best-fit model was found inside a grid prepared by vary- 
ing the parameters in the following intervals: 5 < T(K) < 50 in 
steps of 1 K; < 8 < 5 in steps of 0.5; and 10 < A Q (j^m) < 600 
in steps of 10 fim^ These large intervals, larger than physi- 
cally plausible, were chosen to ensure that the best-fit param- 
eters would not fall on the border of the grid. The best fit was 
then chosen as the one with the smallest^ 2 . During the search, 
we have applied a few constraints: a) T b > r ? , i.e., the grey- 
body envelope must be colder than the blackbody component; 

b) Oft < Q g , the blackbody must be smaller than the greybody; 

c) the model must predict a 24 fim flux smaller than the upper 
limit; d) the predicted flux at 1.1 mm must be smaller than the 
measured value; and e) the envelope mass must be in the range 
of 0.1-10 M . 



2 We took into account the quantization of the parameters in the grid: 
for instance, all the good models for Bl-bS have p = 2.0, and, since the 
step in j} is 0.5, we quote the result as /? = 2.00 ± 0.25. 



where a is the sound speed corresponding to the temperature 
found from the fit, assuming a molecular weight ft = 2.33/ih = 
3.90 x 10~ 24 g. For R we used the radius found from the fit, i.e., 
R = D y/Clg/n. 



A. 1 . Photometry corrections 

Before comparing the model fluxes with the observed data, two 
steps were performed. For the first step, a colour correction was 
made; the flux calibration of PACS and SPIRE was performed 
under the usual assumption that the SED of a source displays a 
flat vF v spectrum. For all other kinds of SED, the derived fluxes 
must be colour corrected according to the intrinsic source spec- 
trum. To correct the fluxes properly, however, we need to know 
the spectrum a priori but that is actually what we want to derive 
from the data themselves. To overcome this circular problem, 
we computed a correction the other way round: a set of grey- 
body models were computed over a large frequency range and 
the fluxes were derived as 



F y (T,^o,^)RSRF v dv 



(A.4) 



where RSRFv stands for the relative spectral response func- 
tion of each PACS or SPIRE filters, computed by the instru- 
ments control teams. To obtain the corresponding flux density, 
F cc was divided by the appropriate filter width that we com- 
puted by imposing the condition that the colour corrections at 
the effective wavelengths is 1 for a SED of constant vF v . The 
derived AA and Av are reported in Table |A.1| As a consistency 
check, we compared the colour corrections we obtained for a 
blackbody with the values found by the instrument teams (see 
http://herschel.esac.esa.int/). For PACS, the agreement is better 
than 1% in all the bands for T > 7 K. The agreement at lower 
temperatures is the same at 100 fim and 160 fim, while it starts 
to diverge at 70 fim, in particular our correction for T — 6 K 
is 1.3% higher, at T — 5 K is 37% higher. At such extremly 
low temperatures, however, the colour corrections are quite dif- 
ficult to be estimated since the intrinsic spectrum differs notably 
from the calibration spectrum. For SPIRE the comparison is less 
obvious because the colour corrections for a blackbody are not 
tabulated as a function of the temperature, but only in the lim- 
iting case of a blackbody in the Rayleigh-Jeans regime. For a 
blackbody at 100 K, we found that our colour correction is less 
than 1.3% higher at 250 fim, and less than 1% for the other two 
bands. 

The colour corrections for the attenuated blackbody were 
found by inverting the definition of the greybody, i.e., 
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Table A.l. Effective wavelengths and filter widths used to com- 
pute the colour corrections. 



A Qim) 


70 


100 


160 


250 


350 


500 


AA (jim) 
Av ( 10" Hz) 


10.6 
6.48 


17.0 
5.09 


30.2 
3.54 


46.3 
2.22 


60.5 
1.93 


113 
1.35 



Table A.2. Colour correction factors for the best fit models re- 
ported in Table[5] Herschel fluxes reported in Fig.[3]are the mea- 
sured fluxes given in Table [2] divided by these factors. 



Source 


70 


100 


160 


250 


350 


500 


Bl-bS 
Bl-bN 


0.790 
0.744 


0.978 
0.946 


0.999 
0.909 


1.016 
1.029 


1.006 
1.015 


0.987 
0.992 



B v (T)e- T = B V (T) - G V (T), from which (B v (r)e- T )cc = B CC (T) - 
G CC (T). In this way, the same grid of models could be used. 

Once the best fit was found, we computed the true colour 
correction factors, i.e., the amount by which the observed fluxes 
must be multiplied to get the instrumental corrected fluxes. For 
the best fit models reported in Table [3] the colour corrections / cc 
are given in Table |A.2| 

For the second step we made before comparing model fluxes 
with observed fluxes, we took into account the Gaussian fit used 
to derive the photometric flux. For example, it is known that the 
PACS PSF is not a Gaussian, so an error is introduced in the 
photometry because of the Gaussian fit. To derive these errors, 
we reduced a set of PACS observations of flux calibrators and 
the results of aperture photometry and of synthetic photometry 
(Gaussian fitting) were compared. This exercise was repeated for 
a set of isolated compact sources in the Perseus field with differ- 
ent integrated fluxes. For 70 /mi, 100 //m, and 160 /mi, we found 
correction factors of 1.6, 1.5, and 1.4, respectively. We are mak- 
ing additional tests to improve our knowledge of Gaussian fit 
errors. For SPIRE, the PSF is much closer to a Gaussian profile, 
so that we do not need correction factors for SPIRE fluxes. 

In summary, Table [2] reports the measured fluxes multiplied 
by the factors that correct for the Gaussian fit. In Fig. [3] the 
Herschel data are the fluxes from Table [2]multiplied by the color 
correction. For instance, the flux at 70 fim of Bl-bS resulting 
from the Gaussian fit is 0.138 Jy, which becomes 0.22 Jy after 
the multiplication by 1.6. This is the value reported in Table [2] 
After being multiplied by 0.79, see Table |A2j the flux becomes 
0.17 Jy, which is the value used in Fig. [3] 

We do not yet have estimates of the uncertainties in the maps, 
so it is not possible to give reliable uncertainties to each flux. The 
uncertainties given in Table [2] were found after a comparison of 
CuTEx fluxes with those found with the getsources algorithm 
( |Men'shchikov et al.|20T2"l >. 

A.2. Photometry in SCUBA maps 

The effective beam of SCUBA at 450 fim is a Gaussian of 
FWHM 17'.' 3 ( [Pi Francesco et al.||2008] l, but the actual beam 
consists of two different spatial components: an inner one with 
a FWHM of 11", and an outer one with a FWHM of 40". The 
sizes we derived for Bl-bS and Bl-bN are smaller than 17'.'3. 
Since CuTEx is more sensitive to compact sources than to ex- 
tended sources, we conclude that what CuTEx fitted was just the 
inner component of the sources, while the extended component, 
larger than the sources' separation, was seen by CuTEx as back- 
ground. To estimate the total flux of the individual sources, we 



used the following approach. If a source has an intrinsic size 6, 
then 



where 6\ = y0 2 + 1 1 2 and similarly for 64. From the analysis 
of Di Francesco et al. (2008 ), we know that the peak fluxes are 
Pi = 0.88P and P 2 = 0.12P, where P is the peak flux of the 17'.' 3 
FWHM Gaussian, from which P 2 = Pi x 0.12/0.88. Finally, 
expressing 84 as a function of 9\ , we have 



, 0.12 

Plot = P\ (01 /ll) + P\ 

iv 1/ j og8 



1 



- II 2 +40 



2\ 



40 2 



(A.5) 



At 850 fim, CuTEx finds sizes that are larger than the effec- 
tive beam of 22'.'9, but since this beam is again the combination 
of two Gaussians ( |Di Francesco et aTT||2008) > of 1975 and 40" 
FWHM, it is likely that also in this case CuTEx fits just the inner 
Gaussian and sees the second one as a background. The formula 
we used to find the total flux is the same as Equation (A.5 I, with 



1975 instead of 11". Since our sources are slightly larger than 
the effective beam, however, the applicability of Equation ( |A.5| l 
to the 850 fim flux is less robust. 

A.3. The mass of a greybody 

Without Herschel data to define the SED in the far-infrared, the 
emitting mass is usually derived from the flux measured at long 
wavelengths where the envelope becomes optically thin. In this 
limit Equation (A.l 1 then becomes 



F v * tB v (T)Q. 



(A.6) 



If the envelope is optically thin, we see the whole mass distribu- 
tion M, so that 



= « J pis - 



(A.l) 



where A is the projected area and K re f is the opacity at the ref- 
erence wavelength A le f. As in other papers of the Gould Belt 
consortium, e.g. , |Konyves et al.] (2010i we adopted an opacity 
of Kmf = 0.1 cmV 1 at ^ ref = 300 jum ( jBeckwith et~aL|l990 L 
With t his choice, the G BS opacity law is very similar to the 



one by Hildebrand ( 1983 1, who adopted Ar re f = 0.1cm g at 



/i re f = 250 fim and B = 2 for A > 250 fim. For this work, actu- 
ally, only K re f is important because B was derived directly from 
the fit. /l re f was shifted to 300 fim for consistency with older re- 
sults obtained with ground-based facilities, like IRAM or JCMT 
By comparing the column density map of IC 5 146 obtained from 
Herschel observations and the GBS dust opacity law with that 
obtained from SCUBA (sub)mm data and with a different opac- 
ity law, Arzoumanian et al. ( 201 1] > found that the GBS dust opac- 
ity law works well for the Herschel data. The exact value for K le f 
is in any case uncertain and likely dependent on the dust temper- 
ature. An uncertainty in the derived mass of a factor ~2 cannot 
be excluded. 



Substituting Equation (A. 7 1 into Equation (A.6 1 and writing 
the solid angle as A/D 2 = nR 2 /D 2 , with R equal to the outer 
radius of the envelope, assumed spherical, and D equal to the 
distance to the source, we arrive at the well known formula 



M ; 



FyD 2 



R,T, (A - 8) 
KyBy(T) 

Since from our best-fit models Aq and B are also found, we de- 
rived a different formula to estimate the mass of a greybody en- 
velope. 



10 



Stefano Pezzuto et al.: Bl-bS and Bl-bN: two first hydrostatic cores candidates in the Perseus star-forming cloud 



At wavelengths where the envelope is optically thin, we find 
by combining Equations (A. 2 1 and ( A.7 1 that 



from which 



M : 



Z? 2 Q / Ao 



Kref Wref 



(A.9) 



This new equation, equivalent for a greybody to the unnumbered 
equation on page 274 in Hildebrand (1983 ), gives the mass in- 



side a sphere of radius R that, for the given opacity law, has opti- 
cal depth 1 at Aq. Of course, Equation ( |A.9[ ) can be applied only 
if we know [3 and Aq, but these are what Herschel data allow us 
to find. The mass of the envelopes, discussed in Section 3, were 
found using the above formula where the greybody parameters 
are those derived from the fit and reported in Table [3] 

Appendix B: Sigma clipping with a threshold 
dependent on the sample size 

Before generating the final maps, the bolometers' time series 
have to be corrected for the presence of glitches caused by the 
impact of high energy particles falling onto the detectors. To 
achieve this correction, we exploited the spatial redundancy pro- 
vided by the telescope movement, which ensures that each sky 
pixel of the final map has been observed several times with 
different bolometers. Outlier detection was then made with a 
sigma-clipping algorithm. Namely, given a sample of N values, 
estimates for the mean and standard deviations are derived. All 
the values that differ from the mean by more than a standard 
deviations are considered outliers and removed from the sample. 

One problem with the sigma-clipping algorithm is the choice 
of the number of standard deviations above which a datum is 
considered an outlier. To avoid making an arbitrary choice of a, 
a formula has been derived that makes a dependent on the size 
of the sample. Below, we give some details of the method. 

For a Gaussian distribution with mean m and standard devi- 
ation cr, the probability to get a value x such that m - acr < x < 
m + acr, is given by erf '(a/ v2) where erf(y) is the error function 



erf(y) 



2 



For example, if a = 3, the probability is erf(3/ V2) = 0.9973. 
Obviously, l-erf(a7 V2) gives the probability to get a value 
outside the same interval. This probability tends to zero when 
a — » oo (by definition erf(+oo)=l). 

In a sample of N data, the probability given by the error func- 
tion, multiplied by N, can be interpreted as a the expected num- 
ber of points inside or outside a certain interval around the mean. 
For any value of a, the number p of points that differ more than 
acr from the mean is then 



p(a; AO 



Hi 



■N 



For an ideal Gaussian, any value of a is allowed. In a real exper- 
iment, however, p is an integer number and data differing from 
the mean by acr, such that p(a;N) <K 1, are suspicious. To im- 
plement this condition in the sigma-clipping algorithm, we have 
to define a precise value of a, say a, above which observed data 




2.0 2.5 3.0 3.5 4.0 

Number of a 



Fig. B.l. In this figure we show, as a function of the number of 
standard deviations, the ratio between N = N{a) as given by 



Equation (B.l I, and N = N(a) approximated with the parabola 
given in Equation (B.2 1. The difference between the two curves 
is always less than 3%. Fitting a cubic instead of a parabola does 
not improve significantly the approximation and makes the rela- 
tion N = N(a) more nuisance to invert. 



are considered outliers and removed from the sample. To this 
aim, we define a as 



p(a;N) = 1 



erf(— ) .1-1 
lV2/ N 



(B.l) 



All the points, if any, that are outside the interval m - acr < 
x < m+acr are considered outliers and removed from the sample. 
This is equivalent to assume that p(a > a;N) = instead of 
p(a > &;N) < I. 

When N ^ oo then a — > oo too. The mathematical property 
of the Gaussian distribution that for an infinitely large sample 
any value is allowed is preserved. 

" for 



It is not possible to invert analytically Equation (B.l 



a given N, but it is possible to derive an approximate analyti- 
cal expression. Indeed, we found that in the range 1 < a < 5 
Equation (B.l i, written in terms of log(A0, can be approximated 



with the parabola 

log(A0 = 0.0794 + 0.2282a + 0.2005a 2 



(B.2) 



Equation (B.l 
than 3%. 



In Fig. B.l we show the ratio between N as given by 



and by Equation (B.2i. The ratio is always less 



Inverting Equation (B.2 1 is trivial, the result is 



-0.569 + V-0.072+4.991og(A0 



(B.3) 



For a given N, the value of a is found from the above equa- 
tion. The number of expected points distant from the mean by 
more than acr is zero. If one or more outliers are instead found, 
they are removed. New values of m and cr are recomputed and 
the new N' < N is used in Equation (B.3 i. This process is re- 
peated until no other outliers are present in the sample. The 
procedure may not converge, for instance, if the noise is not 
Gaussian. For this reason, we do not actually iterate the formula. 
Instead, for each point of the map the detection of outliers is 
done only once. The number of outliers found at the first itera- 
tion, however, can be larger than 1 . 
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In the central region of the map, where the coverage (and 
thus N) is high, the number of standard deviations for clipping is 
larger than in the outskirts of the map where the coverage is low. 
For instance, if a sky pixel has been observed with 40 bolome- 
ters, the above formula gives a = 2.25. So, once we have esti- 
mated the mean m and the standard deviations cr, all of the values 
Xi such that ABS(x, - m) > 2.25cr are flagged as outliers. If in- 
stead a pixel has been observed with 20 bolometers the threshold 
lowers to l.96cr. 
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Fig. 1. Maps of E5xE5 of Bl-bS and Bl-bN at all wavelengths, with sizes and positions derived from the Gaussian fits. The 
FWHM of the respective instrumental beam is shown as a circle in the bottom-left corner of each respective panel. At A < 100 /mi, 
the source S295 is also visible. Colour bars are in Jy/beam for the SCUBA maps, and MJy/sr in all the others. In the Spitzer images, 
positions and sizes of the sources are copied from the 70 yum map. IRAC maps are in the top row, from left to right: a) 3.8 /mi; b) 
4.5 /mi; c) 5.8 /mi; d) 8.0 /mi. Second row: e) Spitzer 24 /mi; f) PACS 70 /mi, Bl-bN is not really detected and its flux corresponds 
to a 1 a rms point source; g) PACS 100 /mi; h) PACS 160 /mi, S295 is no longer visible. Third row: i) SPIRE 250 /mi; j) SPIRE 
350 /mi; k) SCUBA 450 /mi, where the beam shown is 11", see Appendix A.2; 1) SPIRE 500 /mi. Bottom row: m) SCUBA 850 /mi; 
n) Bolocam 1 . 1 mm, where only one source has been fitted. 
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Fig.3. The best-fit two-component models for Bl-bS and Bl-bN. Filled circles are Herschel data (for Bl-bN the 70 /urn flux is an 
upper limit), open circles are SCUBA data, and upper limits at A < 24 yum are from Spitzer IRAC and MIPS. The solid line is the 
best-fit model, a sum of a blackbody (dashed line dominating at shorter wavelengths) plus a greybody (dashed line more important 
at longer wavelengths). Herschel fluxes are those reported in Table [2] multiplied by the colour correction factors, see Appendix A.l. 
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